
Effect of intracellular loop 3 on intrinsic dynamics 
of human p 2 -adrenergic receptor 

Ozcan et al. 



BioMed Central 



Ozcan et al. BMC Structural Biology 201 3, 1 3:29 
http://www.biomedcentral.com/1472-6807/13/29 



Ozcan et al. BMC Structural Biology 2013, 13:29 
http://www.biomedcentral.eom/1 472-6807/1 3/29 



Structural Biology 



RESEARCH ARTICLE Open Access 



Effect of intracellular loop 3 on intrinsic dynamics 
of human p 2 -adrenergic receptor 

Ozer Ozcan 1t , Arzu Uyar 2t , Pemra Doruker 2 and Ebru Demet Akten 3 * 



Abstract 

Background: To understand the effect of the long intracellular loop 3 (ICL3) on the intrinsic dynamics of human 
|3 2 -adrenergic receptor, molecular dynamics (MD) simulations were performed on two different models, both of 
which were based on the inactive crystal structure in complex with carazolol (after removal of carazolol and 
T4-lysozyme). In the so-called loop model, the ICL3 region that is missing in available crystal structures was modeled 
as an unstructured loop of 32-residues length, whereas in the clipped model, the two open ends were covalently 
bonded to each other. The latter model without ICL3 was taken as a reference, which has also been commonly 
used in recent computational studies. Each model was embedded into POPC bilayer membrane with explicit water 
and subjected to a 1 Lis molecular dynamics (MD) simulation at 310 K. 

Results: After around 600 ns, the loop model started a transition to a "very inactive" conformation, which is 
characterized by a further movement of the intracellular half of transmembrane helix 6 (TM6) towards the receptor 
core, and a close packing of ICL3 underneath the membrane completely blocking the G-protein's binding site. 
Concurrently, the binding site at the extracellular part of the receptor expanded slightly with the Ser207-Asp1 13 
distance increasing to 18 A from 1 1 A, which was further elaborated by docking studies. 

Conclusions: The essential dynamics analysis indicated a strong coupling between the extracellular and 
intracellular parts of the intact receptor, implicating a functional relevance for allosteric regulation. In contrast, no 
such transition to the "very inactive" state, nor any structural correlation, was observed in the clipped model without 
ICL3. Furthermore, elastic network analysis using different conformers for the loop model indicated a consistent 
picture on the specific ICL3 conformational change being driven by global modes. 

Keywords: ICL3, Molecular dynamics simulation, Transmembrane helix 6, G-protein binding site, Ligand docking, 
Essential dynamics 



Background 

As the largest family of membrane proteins in the human 
genome, the G protein coupled receptors (GPCRs) are 
structurally characterized by the presence of seven 
membrane-spanning a-helical segments with an extracel- 
lular N terminus and an intracellular C terminus. Upon 
binding to agonists, a series of conformational changes 
propagate along transmembrane helices and reach the 
intracellular part of the receptor, which directly interacts 
with the hetero-trimeric G-protein. Consequently, G pro- 
teins activation triggers different cascades of events 
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depending on the type of agonists bound to the receptor. 
Therefore, as the initiation point to the flow of signals into 
cells, GPCRs are associated with a plenty of diseases that 
make members of this family significant pharmacological 
targets. 

The first solved X-ray crystal structure of GPCR belongs 
to bovine rhodopsin [1,2], which is followed by the crystal 
structure of human |3 2 -adrenergic receptor (|3 2 AR) in the 
inactive state [3,4]. Since 2007, the cholesterol bound form 
of |3 2 AR (PDB:3D4S) [5], the structure of turkey p^adren- 
ergic receptor (PDB:2VT4) [6], the structure of a methyl- 
ated |3 2 AR (PDB:3KJ6) [7] and various forms of inactive 
states of (3 2 AR bound to antagonists such as ICI 118,551 
and alprenolol (PDB:3NY8,3NY9,3NYA,3PDS) [8,9] have 
been reported. Finally, the nanobody-stabilized active state 
of (3 2 AR in complex with G-protein, has been solved by 
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Rasmussen and his coworkers (PDB:3SN6) [10,11]. Still, 
these static pictures of the receptor remain insufficient to 
describe the dynamic character of the receptor, which gov- 
erns the function. It is a well-established concept that pro- 
teins have an intrinsic ability to sample an ensemble of 
distinct conformations in order to perform certain func- 
tions [12]. The ligand simply selects the optimal receptor 
conformation for binding followed by an induced fit to 
stabilize the final conformation. Many questions remain 
on these multiple, ligand-specific conformational states of 
(3 2 AR with different levels of activity from fully active to 
fully inactive, which induce distinct signaling pathways. 

The ternary complex model proposed in 1980 by Lefko- 
witz and his coworkers [13] describes an allosteric mech- 
anism for receptor activation. The agonist molecule, when 
bound to the extracellular part, simply promotes and sta- 
bilizes the high affinity |3 2 AR-G protein complex. Follow- 
ing the laws of thermodynamics, binding of G -protein 
increases the receptor s affinity for agonist binding to the 
same extent. Fluorescence spectroscopic studies of p 2 AR 
by Ghanouni et al [14] presented a model with multiple, 
agonist-specific receptor states, in which the activation oc- 
curs through a sequence of conformational changes. They 
also suggested that the activation barrier for transition 
from intermediate to active state is high, and that in vivo 
the barrier is more likely reduced by G protein binding. 
The presence of an intermediate state is further supported 
by the fluorescence spectroscopy studies of Swaminath 
et al [15,16], suggesting a mechanistic model for GPCR 
activation, where agonist binding stabilizes a series of con- 
formational states with distinct cellular functions. 

In addition to experiments, several MD simulation stud- 
ies have been conducted after the inactive and active states 
of the receptor have been solved by X-ray crystallography. 
One simulation study by Dror et al [17] reveals that the re- 
ceptor exists between two distinct inactive conformations 
of the receptor, one with the ionic lock intact and one with 
the lock broken. In 2011, Dror and his coworkers proposed 
a completely different activation mechanism in which the 
structural changes start at the G protein binding site propa- 
gating upwards as opposed to agonist-induced conform- 
ational changes that start at the agonist binding site and 
propagate down to G protein binding site [18]. The 
agonist-bound crystal structure of |3 2 AR without a binding 
partner (PDB:3PDS) recently revealed by Rosenbaum et al 
[9] is found to be identical to the inactive state of the recep- 
tor (PDB:2RH1). This suggests that in the absence of a 
G-protein, the receptor prefers to adopt the inactive con- 
formation whether or not it is bound to an agonist. In other 
words, the agonist molecule is not sufficient alone to shift 
the equilibrium to the active state. Dror et al [18] also pro- 
posed an intermediate state for G-protein binding site, 
which exists as a part of the receptors intrinsic dynamics. 
Binding of a G-protein to this binding site simply promotes 



a transition to the active conformation, which is further sta- 
bilized by an agonist bound at the extracellular region. The 
most important feature about the dynamics of (3 2 AR is the 
strong coupling that exists between the intracellular G- 
protein binding site and the extracellular ligand-binding site 
of the receptor [7,19]. The receptor behaves like a pair of 
pincers where the intracellular part becomes narrower as 
the extracellular part becomes wider, and vice versa. 

Due to its unstructured nature, ICL3 region is either un- 
resolved in crystallographic experiments or completely re- 
moved and replaced by T4-lysozyme (T4L) to facilitate the 
crystallization. Thus, none of the experimental and simu- 
lation studies have discussed the possible effect of ICL3 on 
the intrinsic dynamics of the receptor. Its replacement by 
T4L to facilitate crystallization did not prevent agonist- 
induced conformational changes based on fluorescence 
spectroscopy measurements [20]. Yet, it is well accepted 
that its direct interaction with G-protein probably have a 
significant role on the receptor s dynamics and the activa- 
tion/inactivation pathway [21,22]. 

In this study, the effect of ICL3 on receptor s conform- 
ational dynamics was investigated via two distinct models 
of the receptor. Both models were generated from the in- 
active state of the receptor (PDB:2RH1) after removal of 
T4L. Moreover, the partial inverse-agonist carazolol was 
removed from the binding site of both models, since the 
goal of this work was to provide data about the intrinsic 
dynamics of the receptor, i.e., the ensemble of conforma- 
tions accessible to its apo form. According to the current 
view on ligand binding, the equilibrium distribution of 
conformational states may be shifted upon ligand binding. 

In the so-called loop model, the ICL3 region was mod- 
eled as an unstructured loop of 32-residues length and 
inserted between two open ends of TM5 and TM6 (resi- 
dues 230 and 263), whereas in the second model, these two 
open ends were "clipped" or simply covalently attached to 
each other. The "clipped" model of the receptor, serving as 
a reference in our study, has been commonly used in recent 
simulation studies as well [18,23]. Both models were sub- 
jected to 1 us MD simulation in a POPC membrane bilayer 
at 310 K. The essential dynamics analysis was carried out to 
reveal important allosteric coupling within the receptor in 
the presence of ICL3. Two distinct snapshots taken from 
the loop models trajectory were further used as docking 
targets for an agonist and an antagonist molecule in order 
to investigate the effect of ICL3 on binding site conforma- 
tions. Finally, elastic network analysis was performed on 
different conformations and loop models to reveal a con- 
sistent picture on receptor intrinsic dynamics. 

Results and discussion 

Presence of ICL3 affects RMSDs and loop mobility 

The difference between loop and clipped model dynam- 
ics is illustrated in the root mean square deviation 



Ozcan et al. BMC Structural Biology 2013, 13:29 
http://www.biomedcentral.eom/1 472-6807/1 3/29 



Page 3 of 16 



(RMSD) profiles throughout the simulation. In Figure la, 
three different RMSD values are plotted for the loop 
model after fitting all snapshots to the initial snap- 
shot based on coordinates of either the whole protein 
(red line: ALL Loop), the core region excluding ICL3 
(green line: CORE Loop) or the transmembrane re- 
gion composed of helices (blue line: TMEMB Loop). The 
large RMSDs observed in loop model (ALL) are due to 
the presence of ICL3. The core region of the loop model 
is equivalent to the clipped model and thus, the com- 
parison between their RMSD profiles (CORE Loop and 
ALL Clipped) reveals that the clipped model has reached 
a plateau at an earlier time (-200 ns) than the loop 
model (-700 ns). The time at which the loop models 
core region reaches the plateau corresponds to the time 
when ICL3 stabilizes as well. The stabilization of ICL3 
corresponds to a significant change in its conformation, 
observed as close packing underneath the receptor, 
which will be discussed later in the text. 

On the other hand, the RMSD values of the trans- 
membrane region (TMEMB) reach a plateau at around 
50 ns in both models. This indicates that the structure 



of transmembrane region is preserved comparably in 
both models. Furthermore, the presence of ICL3 affects 
the mobility of small intra- and extracellular loops 
(ICL1, ICL2, ECL1, ECL2 and ECL3, see Figure 2 
described in Methods) in the loop model, which are 
present in the CORE profile but not the TMEMB. There 
is a strong correlation between the RMSD profiles of the 
whole protein and ICL3 (only), given by the uppermost 
red and green lines in panels A and B, respectively. The 
RMSD value of ICL3 is obtained after alignment of the 
core region to the initial structure (RMSD-ICL3/FIT- 
CORE; green line in Figure lb). Thus, the extremely 
high RMSD of whole receptor with ICL3 is a conse- 
quence the high mobility of the long intracellular loop. 

The root mean square fluctuation profiles (RMSF) are 
plotted using the time range as [50 ns-1000 ns] for both 
models (Figure 3). The RMSF of each alpha-carbon atom 
in the protein is calculated based on the average struc- 
ture of the aligned snapshots. The first 50 ns are ex- 
cluded as the equilibration stage of the transmembrane 
region in both models (see Figure la). The extent of the 
average fluctuation during 950 ns is found to be higher 
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Figure 1 RMSD profiles for the loop and clipped model simulations, (a) RMSDs plotted based on alignments using different regions (ALL, 
CORE: without ICL3, TMEMB: transmembrane helices) of the receptor, (b) RMSD of ICL3 only (green line). 
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(a) Loop model (b) Clipped model 

Figure 2 Two (3 2 AR models used in MD simulations, (a) Loop 
model with the missing ICL3 modelled using MODELLER, (b) 
Clipped model without ICL3. 



in the loop model (red line) than that of the clipped 
model (blue line), in almost all protein regions, except 
the ECL3 region. Another time range for the loop model 
was taken as [700-1000 ns] (green line; LAST 300), 
where the whole protein has reached a plateau (see 
Figure lb). 

The most significant difference between two time 
ranges for the loop model is in the mobility of the ICL3 
region, which fell down to 2-3 A in the second time 
range from 16 A (out of the range of Figure 3). In 
addition, a relatively lower decrease in RMSF is observed 
in all parts of the protein including the loops and the 
more stable helices in the second time range. However, 
such a difference cannot be observed in the mobility of 
the clipped model based on the two time frames (not 
shown). These results indicate that the fluctuation of 
ICL3 region in the loop model is directly reflected on 
every part of the protein structure, including the trans- 
membrane regions. Once ICL3 becomes closely packed 
under the receptor at around 700 ns (see next section), 




Figure 3 RMSF profiles of residues in the loop and clipped 

models. Fluctuations are calculated using the whole 1 us trajectory 

(ALL) or the time frames of 0.7-1 \is (LAST 300). Presence of ICL3 

imparts mobility to other loops of the receptor. Packing of ICL3 

under the receptor hinders fluctuations after 700 ns. 
\ J 



the mobility of the transmembrane region decreases 
slightly and becomes more similar to that of the clipped 
model. 

Another important observation about the RMSF pro- 
files in Figure 3 is the mobility of the ECL2 loop region, 
which is at the extracellular side of the membrane and 
plays an important role as an access point to the binding 
site. For the loop model, during the last most stable 
300 ns, the RMSF of ECL2 decreases to 5 A from 6.5 A, 
as a consequence of the decrease in the mobility of 
ICL3. But still, the mobility of ECL2 in the loop model is 
higher than that in the clipped model irrespective of the 
time ranges considered. The higher ECL2 mobility al- 
lows a wider range of conformational sampling, which 
would include the open/closed forms of the gateway to 
the binding site, making the loop models binding site 
more accessible and accommodating for diffusing li- 
gands than the clipped model. 

The conformational change of ICL3 gives rise to a "very 
inactive" state of the receptor 

Figure 4A shows the RMSD profiles of the sixth trans- 
membrane helix (TM6) from its inactive (PDB:2RH1) 
and its active states (PDB:3SN6) in reported crystal 
structures. The RMSD is calculated for the intracellular 
part of the helix composed of residues 267-282. In the 
loop model (red), the deviation from both inactive and 
active states starts to increase at around 600 ns and 
levels off around 800 ns, amounting to a change of 2 A. 
On the other hand, no significant change is observed in 
the clipped model (green). The deviation is illustrated in 
Figure 4B, where the first snapshot of TM6 (blue) is 
close to the inactive state (green) and the last snapshot 
of TM6 (yellow) is found to be away from both the in- 
active and the wide-open active state (magenta). The 
second view of the receptor from the intracellular side 
in Figure 4B shows that ICL3 becomes more wrapped 
up under the core of the receptor (final frame shown), 
which will be named as a "very inactive" state inaccess- 
ible to G-protein binding. 

This close packing of ICL3 is quantitatively repre- 
sented in Figure 4C, which shows the sudden change in 
the ^-offset with respect to the j-offset. The value of x- 
offset is the difference in the x coordinates between the 
center of masses for the core of the receptor and the 
ICL3 region. Similarly, the value of j-offset is calculated 
from the difference in y coordinates of these center of 
masses. Three distinct clusters are observed in time 
ranges of [0-470] ns, [470-670] ns and [670-1000] ns. 
The first (red) and third clusters (blue) correspond to 
the open and packed states of ICL3, respectively. The 
second one (green) represents a transition between the 
two states. Interestingly, the second clusters starting 
time at 470 ns corresponds approximately to the time at 
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Figure 4 Conformational change in the intracellular part of TM6. (A) RMSD profile of the TM6 intracellular part from its inactive and active 
states in the loop (red) and clipped (green) models. (B) Initial (blue) and final (yellow) snapshots from 1 us loop trajectory, showing TM5 and 
TM6 in comparison with the inactive (green) and active (magenta) crystal structures. (C) Changes in the x-offset with respect to the y-offset 

representing the change between center of masses of ICL3 and protein's core. 

^ J 



which small helical formations appear inside ICL3 (see 
Additional file 1: Figure SI). It is an open question 
whether these helical formations might trigger the tran- 
sition to the packed state. Additionally, the onset of the 
third cluster around 670 ns corresponds to the time, at 
which TM6 starts to deviate from the reference crystal 
structures (see Figure 4A). 

The changes observed at the intracellular part of the 
receptor seem to affect the extracellular part, specifically 
the binding site of the receptor. Figure 5a illustrates the 
profile of the distance between the pair of residues, 
Ser207-OG on TM5 and Aspll3-CG on TM3. Recent 
simulation studies by Katritch et al. [24] have revealed 
that tilting of TM5 towards the receptor axis enables an 
optimum interaction between agonists and the two an- 
chor sites, Aspll3/Asn312 and Ser203/Ser204/Ser207 
side chains. Based on experimental studies [25], the dis- 
tance between the side chain oxygen of Ser207 and 
gamma carbon of Asp 113 should be within a range of 
8 A (blue, horizontal line) and 10 A (purple line) in 
order to accommodate the agonists at the binding site. 
However, the distance profile in Figure 5a is most often 
out of this critical range and even beyond the distance 



values of the inactive state, which is around 11 A. At 
around 600 ns, when the sudden conformational changes 
in both TM6 and ICL3 occur, this distance starts to in- 
crease from 13 A to 16 A in loop model (red). On the other 
hand, there is no significant change in the distance profile 
of the clipped model (green) after 600 ns as expected. 

The variation in the distance between Ser207 and 
Aspll3 is illustrated in Figure 5b (top view, looking 
down from the extracellular side), which shows two con- 
formations of the receptor, with minimum and max- 
imum distance values of 8.3 A and 18.6 A (shown as 
light blue and purple, respectively). Clearly, the intracel- 
lular part of TM5 is slightly moving into the core region 
(see Figure 5b), while its extracellular part is moving 
away from the core region of the receptor. As a result, 
Ser207, which is located at the extracellular part of 
TM5, drifts away from Aspll3 on TM3, position of 
which does not change notably. Similar motions are ob- 
served for TM4 and TM6 as well. As a result, the en- 
larged binding site becomes unfavorable for agonist 
binding due to lack of some key interactions. 

Figure 5c illustrates the correlation between the TM6 
shift at the intracellular part of the receptor and the 
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Figure 5 Conformational change in the binding site, (a) Ser207(Oy)-Asp1 1 3 (C y) distance profiles for the loop (red) and clipped (green) 
models. Blue and purple lines bracket the experimentally determined critical range in active states [25]. (b) Conformations representing the open 
and closed forms of the binding site in loop model, (c, d) The correlations between TM6 shift and Ser207-Asp1 13 distance for the loop and 
clipped models, respectively. 



Ser207-Aspll3 distance at the extracellular part of the 
receptor. In the loop model, the change in the RMSD 
value of TM6 with respect to the active state happens at 
around the same time as the increase in the Ser207- 
Aspll3 distance, whereas in the clipped model, no such 
correlation is observed (see Figure 5d). Three different 
states of the structure are observed in the loop model at 
around [0-640], [640-770] and [770-1000] ns intervals. 
In contrast, the conformational variations of the clipped 
model remain in a restricted area, which corresponds to 
the first conformational state ([0-660] ns) of the loop 
model. 

In addition to the 1 [is MD simulation, three inde- 
pendent 100 ns MD simulations with different initial 
conformations and velocities were performed as ex- 
plained in Methods. Due to restrained conditions during 
the preparation stage prior to MD runs, the Ser203- 
Aspll3 distance value was extended to -16 A from the 
initial value of -11-12 A. Ser203 is another key residue 
for binding that lies on the next turn above Ser207 on 
TM5. The change in the distance between Ser203 and 
Aspll3 is mainly the result of a change in the position 



of Ser203 (both backbone displacement and side chain 
rotation). Within the first 20 ns of all four MD simula- 
tions (including the l\is simulation), the distance rapidly 
decreases back to its initial value of -11-12 A as shown 
in Figure 6. The explanation for such a decrease is that 
all four simulations started with an expanded extracellu- 
lar (binding site region) part and an intracellular part 




0 10 20 30 40 50 60 70 80 90 TOO 



Time (ns) 

Figure 6 Ser203(OY)-Asp1 13(Cy) distance profiles. The initial 

distance indicative of a relatively open binding site decreases in all 

simulations (1 us and three 100 ns independent trajectories) 

resulting from the positioning of TM6 and ICL3. 
v> J 
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with ICL3 set aside. However, the receptor is found to 
be in equilibrium when its extracellular part is open 
(wide) and its intracellular part is closed (narrow) with 
ICL3 closely packed underneath or vice versa. This is 
the direct consequence of the strong allosteric coupling 
that exists between extracellular and intracellular regions 
of the receptor and this seems to be consistent with the 
ternary complex model suggested by de Lean et al [12]. 

The loop model shows a conformational variation in its 
second intracellular loop (ICL2), which is correlated with 
the motion of ICL3. Additional file 2: Figure S2 shows four 
different stages of conformational variation of ICL2 
changes with respect to intracellular part of TM6. Both 
RMSD values are calculated with reference to the active 
state (PDB:3SN6). However, no such major structural 
change is observed in the clipped model. The motion of 
ICL2 in the loop model is also illustrated, in which the 
ICL2 between TM3 and TM4 steps aside as the ICL3 
comes closer to the middle region (initial stage: blue, sec- 
ond stage (at 700 ns): cyan, final stage: red, and active 
crystal: green). No such conformational rearrangement in 
ICL2 is observed in the clipped model. Additional file 3: 
Figure S3 illustrates the change in the RMSD value of 
ECL2 with the change in the Ser207-Aspll3 distance. As 
the distance increases in the loop model, there is a con- 
formational variation in ECL2 with respect to the active 
state. However, no such correlation is observed in the 
clipped model. Clearly, the structural variation in ECL2 is 
directly affected by the change in the distance as a result 
of a shift of TM5 away from the binding site, which is in 
turn a consequence of the ICL3 motion and TM6 shift at 
the intracellular part of the receptor. 

Ionic lock (Arg131-Glu268) is not a molecular switch 

During 1 us long MD simulation of the loop and clipped 
models, which represent the inactive state of the recep- 
tor, the ionic lock profiles are monitored as shown in 
Additional file 4: Figure S4. Consistent with previous work 
[18], the ionic lock seems to be on and off during the sim- 
ulations of both models (in upper panels). Thus, this ionic 
lock cannot distinguish between active and inactive states. 
The cause behind the breakage/formation of this ionic 
lock is found to be the result of a change in the rotational 
state of the X angle of Glu268, which coordinates perfectly 
well with the ionic distance profiles (in lower panels). 

Furthermore, the increase in the distance between two 
side chains that form the ionic lock, namely Argl31-N 
and Glu268-0, coincides properly with the increase in the 
distance between their alpha- carbons. In the profiles of 
the inactive state, the backbone distance fluctuates at 
around 9.5 A and reaches 12.4 A at most. However, in the 
known crystal structure of the active receptor (PDB:3SN6), 
the distance between alpha-carbons is around 16 A as a 
result of a significant outward shift in the intracellular part 



of TM6. Thus this backbone distance could be one pos- 
sible measurement for detection of activation. 

Essential dynamics analysis reveals the transition to the 
"very inactive" state in the first principal mode 

For both models, each frame in the trajectory was 
aligned onto the initial structure. Then principal compo- 
nent analysis (PCA) based on only Ca atom coordinates 
was performed to understand the effect of ICL3 on the 
essential dynamics of the receptor [26]. The first princi- 
pal mode explains about 69% and 22% of the proteins 
overall motion in the loop and clipped models, respect- 
ively. Figure 7a displays three different RMSD profiles 
for the intracellular part of TM6 in the loop model after 
alignment on the active state. These are computed for 
the original trajectory (blue, same profile as in Figure 4A) 
and for the two reconstructed trajectories, one including 
only the first mode (red) and the other showing the cu- 
mulative effect of the first five modes (green). Projection 
of the MD trajectory onto the first principal eigenvector 
(red) shows an abrupt change, which is coupled with the 
transition to the "very inactive" conformation. The pro- 
file obtained from the projection of cumulative five 
modes, as expected, explains this transition better. 

Additionally, the distance profiles between Ser207 and 
Aspll3 in the loop model was recalculated using Ca 
atoms only as shown in Figure 7b. A high correspond- 
ence between the original and the two reconstructed 
profiles is observed as in Figure 7a. Also, a plot of 
RMSD value versus the Ser207-Aspll3 distance clearly 
shows that the essential modes (first and cumulative 
five) describe the distribution in the original trajectory 
(Figure 7c). Thus, the closure of the ICL3 driven by the 
first mode (Figure 7d) is strongly coupled with the 
opening of the binding site indicated by the Ser207- 
Aspll3 distance. For the clipped model, same profiles are 
plotted in Additional file 5: Figure S5. The profiles ob- 
tained from the projection of the first and cumulative five 
modes do not explain satisfactorily the dynamics of the 
extracellular and intracellular regions of the receptor. 

Elastic network modeling reveals coupling between 
global modes and ICL3 conformational transition 

MD simulations were performed on a receptor model 
including a specific unstructured conformation of ICL3 
obtained from MODELLER. In order to show the inde- 
pendence of the observed phenomena from the initial 
MD structure, ANM was performed on four distinct 
conformations of the loop model. These were selected as 
the initial, average and final structures of the 1 us -long 
loop trajectory, and a receptor model containing an al- 
ternative unstructured conformation of ICL3, also pro- 
vided by MODELLER. The RMSD between the 
alternative loop model and the one used in our MD 
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Figure 7 Essential dynamics of the loop model, (a) RMSD profiles for the intracellular part of TM6 in the loop model's original trajectory and 
after the projection onto the first and the cumulative five principal modes, (b) Original and reconstructed profiles for the distance between 
Ser207(Ca) and Aspl 13(Ca). (c) The correlation plot between RMSD of TM6 and Ser207-Asp1 13 distance, (d) Projection of the loop model's 
trajectory onto the first principal mode, shown as harmonic motion. 



simulations is ~ 20 A for the ICL3 residues after an 
alignment of transmembrane regions. 

The correlation between PCA and ANM modes is 
routinely assessed by the average overlap value 



O a 



1 k k 

EE 

, <=i ;=i 



1/2 



(1) 



where and u ; represent the i th and ; th normalized ei- 
genvectors from PCA and ANM, respectively. The 



squared inner dot products are generally summed over the 
first k =10 modes, which describe the collective subspaces of 
each method. The average overlap values are 0.64, 0.72, 0.62 
and 0.66 between the first 10 modes of PCA and ANM per- 
formed with the initial, average and final structures of MD 
run and the alternative loop model, respectively. These 
values are quite high (relatively closer to 1), representing sat- 
isfactory overlap between ANM and PCA subspaces. 

In Figure 8, the conformational changes of ICL3 and 
TM6 are shown for specific ANM modes that yield a 
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Figure 8 Collective modes driving ICL3 and TM6 towards the core of the receptor. (A) The first anharmonic principal mode for the 1 us 
MD trajectory (same red and blue conformers as in 7D, cyan representing the average MD structure), (B) the second slowest mode from ANM 
performed on the average structure, and (C) the first mode from ANM performed on the alternative loop model. On the left panels, alternative 
conformations (red and dark blue) are provided together with the average structure (cyan) for the specific mode. The right panels provide vector 
representations of the same deformations based on average structures (intracellular view). 
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high overlap value (Oy= |pi.U/|) with the first principal 
component of MD run. The left and right panels on the 
figure indicate different views (side and intracellular) for 
each mode. Specifically, vector representations of deforma- 
tions, shown from the intracellular part (right panels), indi- 
cate clearly the inward movement of ICL3 and TM6 in the 
1st mode of MD (Figure 8A), ANM 2nd mode for the aver- 
age structure (0 1>2 = 0.81, Figure 8B) and ANM 1st mode 
for the alternative loop model (0 1?1 = 0.57, Figure 8C). 

The overlap matrices calculated based on residue 
displacements of ICL3 and TM6 region only (see 
Additional file 6: Figure S6) indicate several ANM modes 
exhibiting high overlap with the first mode of MD. Thus, 
slow modes of ANM clearly drive the significant con- 
formational change of ICL3 and TM6 towards the recep- 
tor core, independent of the ICL3 conformation/model 
used. In summary, our ANM analysis justifies that the 
ICL3 dynamics observed in MD run can be attributed to 
be a feature of intrinsic receptor dynamics in conformity 
with a recent study carried on catalytic loop motions for 
different enzymes [27]. 

The clustering of MD snapshots reveals more 
conformational variations in the loop model 

The snapshots taken from the simulations of the loop 
and clipped models are clustered all together based on 



different regions of the receptor: transmembrane region, 
intracellular part of the receptor, ICL2 and ECL2 loop re- 
gions using an RMSD cutoff of 1.8 A, 1.8 A, 3.7 A and 
3.3 A, respectively. The region for the alignment is chosen 
as the transmembrane region in all four cases. In all cluster 
profiles shown in Figure 9, the simulation time is divided 
into five ranges. Frames 1-5000, 5001-10000 and 10001- 
11500 are taken from the (is-long loop run, (is-long clipped 
run and the three short runs for the loop model. 

The clustering profile of the transmembrane region 
shows four distinct clusters for 1 [is simulation of the loop 
model as illustrated in Figure 9a. Two of those clusters 
dominate over the other two, since they contain 54% and 
37% of the total snapshots, which are observed at around 
[0-600] ns and [600-1000] ns, respectively. On the other 
hand, two distinct clusters are obtained for the clipped 
model, and only one of them dominates for 92% of the 
time. For each of the three short MD simulations, there is 
only one single cluster that dominates during 100 ns. 
These results indicate that the transmembrane region of 
the loop model alternates between two distinct conforma- 
tions, while the clipped models transmembrane region 
prefers to adopt only one. Interestingly, the second con- 
formation in the loop model observed between 600 and 
1000 ns coincides with the time at which the ICL3 
changes its conformation and the receptor adopts a "very 
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Figure 9 Clustering profiles of all trajectories based on different regions of the receptor, (a) transmembrane region, (b) G-protein binding 
site, (c) ICL2 and (d) ECL2. All graphical representations were prepared using MATLAB [28] in this article. 
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inactive" state with an expanded binding site. A highly 
similar clustering profile is obtained for the binding site 
region as shown in Additional file 7: Figure S7. This is an 
expected outcome considering that the binding site is em- 
bedded in the transmembrane region. 

Figure 9b shows the cluster profile of the intracellular 
part of the receptor, which consists of residues interacting 
with the G -protein based on the active crystal structure 
(PDB:3SN6) [29]. Clearly, the loop models intracellular 
part samples three distinct states while the clipped models 
intracellular part only samples one conformation. In the 
three short simulations, there is also one single conform- 
ation dominating the others. The structural flexibility of 
the intracellular part is critical in making contact with the 
G protein. For the loop model, three snapshots were se- 
lected from each cluster shown in Figure 9b as representa- 
tives and illustrated in Additional file 8: Figure S8 with a 
bottom view to show the contact of the receptor with the 
helical segment of gamma subunit of G protein. In the ac- 
tive crystal structure (PDB:3SN6) taken as a reference and 
placed on top of the figure, G proteins helical segment 
nicely fits the binding cavity. At the initial stages of the 
simulation, the binding cavity is almost preserved. Towards 
the end of the simulation, the motion of ICL3 closes down 
the G-protein binding site almost completely, as shown in 
the last frame leaving no contact point for the G protein. 

The clustering profile in Figure 9c, shows two dominant 
clusters for the ICL2 region in both model. The two dis- 
tinct states in the loop model are sampled for about 29% 
and 58% of the time, while in the clipped model two major 
clusters are sampled for 39% and 59% of the time. Each of 
the short simulations of the loop model still does not show 
structural variation during 100 ns, similar to other two 
cases above. These results indicate that ICL3 has no sig- 
nificant effect on the conformational sampling of ICL2. 

Finally, the cluster profile in Figure 9d shows four dis- 
tinct clusters for ECL2 region of the loop model, with 
each consisting of a considerable amount of snapshots 
and sampled consecutively during the simulation. On 
the other hand, the clipped models one us simulation as 
well as three short simulations of the loop model, impart 
no conformational variation to the ECL2 region. The 
ECL2 loop region is the second extracellular loop cover- 
ing the top of the receptor and plays a critical role of 
providing a passage to the binding site region. Therefore, 
the ability of ECL2 to sample various conformations, be- 
ing a functionally important feature for the receptor, is 
clearly enhanced in the presence of ICL3. 

Docking results of epinephrine and ICI to an open and a 
closed form of the binding site 

Two frames of the loop model are selected from the 1 
us trajectory to represent the two extreme cases of 
Ser207-Aspll3 distance value (see Figure 5a, b). One of 



the conformers is a closed form with a distance value of 
8.31 A, which is in the range of active states (8-10 A) 
[26], while the other conformer represents an open form 
with the maximum distance value of 18.63 A. The 
docked ligands are a natural agonist epinephrine, and 
ICI, which is an antagonist with a known crystal con- 
formation (PDB:3NY8). The epinephrine is chosen due 
to its relatively small size, and ICI is selected because it 
is a large antagonist with an experimentally determined 
conformation. 

Figure 10a illustrates the poses of epinephrine with 
highest scores docked to open and closed forms of the 
receptor (a top view looking down from the extracellular 
region). Epinephrine is shown as sticks while the key 
residues are in ball-and-sticks representation. The epi- 
nephrines highest score conformation docked to closed 
form (light blue) has more favorable interactions with 
neighboring residues than the highest score conform- 
ation docked to open form (magenta) of the receptor 
(see Figure 10b for an alternative side view to the struc- 
tures). There exist a total of eight neighboring residues, 
which interact with epinephrine in closed form within a 
radius of 3.5 A, namely Aspll3, Vail 14, Ala200, Ser204, 
Ser207, Phe289, Phe290, and Asn293. Seven of them, ex- 
cluding Ala 200, are known to be key interacting resi- 
dues in agonist binding [30] (Additional file 9: see Figure 
S9A, B showing interactions of epinephrine obtained 
from MOE tool) [31]. On the other hand, the best pose 
of epinephrine in the open form is found slightly out of 
the binding site region making interactions with ten resi- 
dues with a distance of less than 3.5 A (Aspll3, Vail 17, 
Phe282, Cys285, Trp286, Phe289, Leu311, Asn312, 
Gly315, Asn318), among which only three (Asp 11 3, 
Phe289 and Asn312) are key residues (see Figure 10a). 
The other seven are not reported as being significant in 
agonist binding. Considering the small size of epineph- 
rine, the closed form is found to be more favorable than 
the open form. 

The docking results of ICI indicate a complimentary 
situation. Due to its large size, the antagonist ICI cannot 
fit into a narrow binding site in the closed form, but it can 
be favorably docked into a wider binding site, such as the 
open form. Figure 10C shows the highest score conforma- 
tions docked in the open (magenta) and closed (light blue) 
form of the receptor (see Figure lOd for a comparison to 
the bound state of ICI in the crystal structure shown as 
black sticks). The RMSD values of docked ICI to open 
and closed forms with respect to its native state are deter- 
mined as 3.95 A and 8.16 A, respectively. Clearly, ICI 
when bound into a wider binding site is able to interact 
with experimentally reported key residues for antagonist 
binding. ICI interacts with eight residues within a radius 
of 3.5 A, Aspll3, Tyrl99, Ser203, Ser204, Phe208, Trp286, 
Tyr308 and Asn312, of which five (Aspll3, Ser203, 
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crystal form 



(c) ICI (top view) (d) ICI (side view) 




Figure 10 Best docked poses of ligands. (a, b) for agonist epinephrine and (c, d) for antagonist ICI in open (magenta) and closed (light blue) 
forms of the receptor, showing both top and side views. Ligands are shown in sticks representation in the same color with its docked structure. 
Key residues interacting with the ligand are labeled and also shown as ball-and-stick. For comparison, the bound state of ICI in the crystal 
structure (PDB:3NY8) is shown in black. All molecular graphics were prepared using PyMOL [32] and VMD [33] in this article. 



Ser204, Trp286, Asn312) are reported as key residues [8] 
(see Additional file 9: Figure S9 for the specific interac- 
tions of ICI). Furthermore, ICI is correctly oriented in the 
binding pocket with the hydrophobic catechol ring inter- 
acting with Ser203 and Ser204, and its polar end interact- 
ing with Aspll3. However, in the closed form, ICI is 
improperly positioned in the binding pocket, interacting 
within a radius of 3.5 A with ThrllO, Aspll3, Vail 14, 
Phel93, Tyrl99, Phe289, Asn293, Lys305 and Tyr308 of 
which only two (Aspll3 and Vail 14) are known key 
residues. 

To sum up, the agonist is favorably bound to the 
closed form, which coincides with the ICL3-open con- 
formation during the first half of the simulation. The 
ICL3-open conformation may correspond to an inter- 
mediate state that promotes G-protein binding, which 
seems to be stabilized by the presence of the agonist at 
the binding site. An MD simulation, where the binding 
site would be constrained to the agonist bound 



geometry, would give more insight in this perspective. In 
the second half of the simulation, the open geometry of 
the binding site is observed to which the antagonist fa- 
vorably docks. This "very inactive" state corresponds to 
the closely packed ICL3 that completely blocks the G- 
protein binding site. 

Conclusions 

The crystal structure of |3 2 AR has been resolved for the 
first time in 2007, and since then in silico studies have 
been conducted to unravel structure-dynamics-function 
relationship of this G-protein coupled receptor. How- 
ever, the missing intracellular loop ICL3, which is known 
to interact with the G protein at the cytoplasmic side, 
has not been considered or elaborated so far in these 
studies. In this respect, our MD study exposed the 
marked effect of ICL3 on collective dynamics and justi- 
fied the correlated motion between the intracellular G- 
protein binding site and the extracellular ligand-binding 
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site of the receptor [7,19]. The receptor behaves like a pair 
of pincers where the intracellular part becomes narrower 
as the extracellular part becomes wider, and vice versa. 

In accordance with this coupling behavior, our us long 
MD simulation of |3 2 AR, which included the modeled 
intracellular loop ICL3, revealed a so-called "very in- 
active" state of the receptor, which has not been reported 
before. In the second half of the simulation, the ICL3 
moved toward the core of the receptor and completely 
blocked the G -protein binding site. Consequently, the 
intracellular part of TM6, adjacent to ICL3, also shifted 
toward the core of the receptor. This conformational 
change in TM6 is in the opposite direction with respect 
to experimentally observed deformation during activa- 
tion [9], which is observed as an expansion or outward 
movement towards the lipid membrane. This new in- 
active state of the receptor may provide insight into the 
design of novel therapeutic drugs. 

Another important observation is the correlated mo- 
tion between the binding site and the G-protein binding 
site regions. At around the same time when the ICL3 
blocked the G-protein binding site, the extracellular 
binding site of the receptor expanded. The expansion 
was detected based on the distance profile between two 
anchor residues, Ser207 and Aspll3, critical for agonist 
binding. In line with these findings, our docking studies 
indicated favorable antagonist binding to the expanded 
binding site (closed ICL3) and agonist binding to the 
closed binding site (open ICL3). 

This coupled conformational change seems to be trans- 
mitted from the intracellular part to the extracellular part 
of the receptor via TM5 and TM6. As Sara Linse stated in 
her review [34], "a 7TM receptor is like a bundle of rods 
immersed in the membrane and if a ligand grips the bun- 
dle at one end, the bundle opens up like a bouquet of roses 
at the other end". In our study, the intracellular end is 
gripped or held tight by ICL3, and consequently the extra- 
cellular part opened up. In other words, ICL3 played the 
dominant role in inducing the change in the intracellular 
part, which induced in turn the extracellular part. This 
dominancy of ICL3 is expected due to its high mobility, 
which is also a desirable quality for initiation of interac- 
tions with intracellular proteins [20,21], Thus, we suggest 
that when left without ICL3, the receptor would not be 
able to sample that inactive state at all. 

This transition to the "very inactive" state took place 
within a time frame of about 0.1 us (starting at ~ 0.6 us). 
In the last 0.3 us of the simulation, ICL3, which was ob- 
served to be the most mobile region of the receptor dur- 
ing the whole simulation, preserved its close state. 
Principal component analysis of 1 us long MD trajectory 
showed that the first principal mode, which explains 
69% of the overall motion governs the transition from 
the initial inactive state to the "very inactive" state. 



At this point it may be argued that our simulation 
conditions, such as the absence of any ligand at the 
binding site and/or intracellular proteins that may inter- 
act with ICL3, depict a non-physiological environment. 
Even though, the receptor function is clearly linked with 
its interacting partners, our aim was to elucidate the in- 
trinsic conformational dynamics of the intact receptor. 
Based on the widely accepted population shift mechan- 
ism [35-37], we tried to uncover the pre-existing con- 
formational states of the apo receptor, which may be 
shifted and/or modified by the presence of binding part- 
ner (s). In fact, our ANM analysis using different con- 
formers/models of ICL3 strengthened our MD results 
on receptor dynamics. Either the first or second collect- 
ive mode in ANM was found to be coupled with the 
specific motion of ICL3, independent of the model used. 
In contrast, none of these conformational transitions, 
nor any allosteric coupling between intra- and extracel- 
lular parts, were observed in the clipped model simula- 
tion lacking ICL3 region. Thus, we stress that the 
presence of ICL3 provides a more realistic constriction 
than those of clipped and non-clipped (loose ends) 
models so far used in |3 2 AR simulations. 

Future works will be focusing on the loop model, 
which will consist of constraining the binding site region 
to observe the reverse transition (or release) from the 
"very inactive" state to the inactive/ intermediate state 
and possibly the active state with an opening of the G- 
protein binding site. 

Methods 

Preparation of the receptor models 

The X-ray crystallographic structure of human (3 2 AR in 
complex with T4 lysozyme (T4L) (PDB:2RH1) at 2.40 A 
resolution [3] was used as the initial conformation. After re- 
moval of T4L, the missing intracellular loop region (ICL3) 
between residues 231 and 262 was modeled as an unstruc- 
tured loop of 32-residues length via MODELLER [38]. Even 
though ICL3 possibly exists as an unstructured loop [33], 
our modeled loop in Figure 2a interestingly resembles a 
conformer generated from a fully extended chain by Dror 
et al [18]. Moreover, it is expected that ICL3 can sample 
various conformations during our 1 us long run. 

In the second model, Leu230 and Lys263, which are 
the two ends of the missing region, were covalently at- 
tached to each other via a peptide bond to form the 
"clipped" model (Figure 2b), which has been commonly 
used in simulation studies. 

Preparation of the environment with the receptor 

Each model was later embedded into a palmitoyloleoyl- 
phosphatidylcholine (POPC) membrane bilayer along 
the z-axis using VMDs Membrane Plug-in Tool, vl.l 
[39]. The receptor was positioned with an oblique angle 
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of 8° between its main principal component along the 
membrane and the z-axis [40]. A total of fifteen internal 
water molecules detected experimentally in the crystal 
structure were retained because they make hydrogen 
bonds with the most conserved residues of the receptor 
and thus possibly contribute to its structural stability. 
Using VMDs solvate module, the protein-lipid system was 
solvated in both intracellular and extracellular sides with a 
thickness of 15 A and 13 A for the "loop" and "clipped" 
model, respectively. Finally, the protein-lipid-water system 
was ionized with Na + and CI" ions to make the total 
charge of the system to be equal to zero, which is ne- 
cessary for Particle-Mesh Ewald summation method used 
in electrostatic energy calculations. The resulting perio- 
dic box dimensions were (86x86x100) and (77x69x90) 
in Angstrom for the "loop" and "clipped" models, 
respectively. 

Molecular dynamics simulations 

Using the suggested procedure for membrane protein sys- 
tem preparation [41], both models were subjected to three 
preparation stages. The system consisted of three compo- 
nents of different types, each having a different response 
time to outside forces. Thus, to reach the equilibrium fast, 
it was practical to keep some components fixed, while 
other components were free to move. The first equilib- 
rium stage consisted of melting the lipid tails where only 
the lipids were free to relax while the protein and waters 
were held fixed. At the end of the simulation, the unrealis- 
tically aligned lipid molecules, transformed into a more 
disordered, liquid-like structure. In the second preparation 
step, the protein s motion was constrained while lipid and 
water molecules were free to move. Finally in the third 
stage, the protein was released, and all components were 
allowed to relax. In every preparation stage, the system 
was subjected to 1000 steps of energy minimization 
followed by 0.5 ns MD simulation. At the end of the third 
stage, the area per lipid was stabilized at 0.635 nm 2 /lipid 
in agreement with the experimentally measured value of 
0.65 nm 2 /lipid [42]. Also, the surface area of the mem- 
brane in -xy directions decreased due to close packing of 
lipid molecules with the protein. 

Each model was later subjected to 1 us MD simulation 
with NAMD v2.7 software tool [41]. CHARMM22 
[43,44] and CHARMM27 [45,46] forcefields were used 
to describe the interaction potential of the protein and 
the lipid respectively, and waters were treated explicitly 
using TIP3P model [47]. The system was composed of a 
total of 68,001 and 42,701 atoms for the "loop" and 
"clipped" model, respectively. MD simulation was per- 
formed at constant NPT at 310 K using Langevin dy- 
namics for all non-hydrogen atoms, with a Langevin 
damping coefficient of 5 ps 1 . The system was kept at a 
constant pressure of 1 atm by using a Nose-Hoover 



Langevin piston [48] with a period of 100 fs and damp- 
ing timescale of 50 ps. Long-range electrostatic interac- 
tions were treated by particle mesh Ewald (PME) 
method, with a grid point density of over 1 A. A cutoff 
of 12 A was used for van der Waals and short-range 
electrostatics interactions with a switching function. 
Time step was set to 2 fs by using SHAKE algorithm for 
bonds involving hydrogens [49] and the data was re- 
corded at every 200 ps. The number of time steps be- 
tween each full electrostatics evaluation was set to 2. 
Short-range non-bonded interactions were calculated at 
every time step. 

For the "loop model", three additional 100 ns long MD 
simulations starting with different initial velocities were 
performed alongside 1 us MD simulation. The aim was 
to possibly explore different conformational subspaces 
than those visited during the long trajectory of the loop 
model. 

Docking calculations 

Docking was performed using the software tool AutoDock 
v4.0 [50]. The docking site was selected based on the loca- 
tion of the partial inverse agonist carazolol in the complex 
structure (PDB:2RH1). Two distinct snapshots taken from 
MD trajectory were used as target structures. Lamarckian 
genetic algorithm was used to explore the conformational 
space. A total of 100 runs were performed for each struc- 
ture with each run consisting of l.OxlO 6 and 1.5xl0 6 
energy evaluations for epinephrine and ICI ligands, re- 
spectively. Grid box constructed with a spacing of 0.375 A 
had dimensions of 24 A x 24 A x 24 A for all dockings. 
For each docking experiment, the pose with the highest 
score (lowest binding energy of AutoDock) was used as 
the most probable solution for that complex. 

Elastic network analysis 

The collective/global modes of the protein were extracted 
via the anisotropic network model (ANM) [51,52], which 
describes the protein as a coarse-grained elastic network 
of harmonic springs based on a minimum-energy folded 
conformation. The network is formed simply, by conn- 
ecting the close-neighboring alpha-carbon atom (called 
nodes) pairs in the folded structure. The slow or the low- 
frequency modes extracted from normal mode analysis of 
the elastic network are known to successfully describe the 
functional conformational changes. 

In our current work, the receptors loop model was 
embedded into a coarse-grained membrane environment 
according to the methodology provided in Lezon et al. 
[53] The membrane consisting of spheres arranged in an 
FCC lattice had a diameter of 80 A and a thickness of 
33 A. The cutoff value for pairwise interactions between 
nodes was taken as 11 A. The force constants of har- 
monic springs were selected as 1.0, 2.0 and 4.0 for 
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protein-protein, protein-membrane and membrane- 
membrane type of pairwise interactions, respectively. In 
this model, the membrane environment serves as a con- 
striction and thereby inhibits the unrealistically large 
fluctuations of the transmembrane helices that would be 
observed if ANM were applied to the protein alone [53]. 
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